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Abstract 



As a mathematical theory for the stochasstic, nonlinear dynamics of in- 
dividuals within a population, Delbriick-Gillespie process (DGP) n{t) G 
Z^, is a birth-death system with state-dependent rates which contain the 
system size F as a natural parameter. For large V, it is intimately re- 
lated to an autonomous, nonlinear ordinary differential equation as well 
as a diffusion process. For nonlinear dynamical systems with multiple at- 
tractors, the quasi-stationary and stationary behavior of such a birth-death 
process can be underestood in terms of a separation of time scales by a 
T* ~ e"^ (a > 0): a relatively fast, intra-basin diffusion for t < T* 
and a much slower inter-basin Markov jump process for t ^ T*. In the 
present paper for one-dimensional systems, we study both stationary be- 
havior (t = oo) in terms of invariant distribution pf^^{V), and finite time 
dynamics in terms of the mean first passsage time (MFPT) T„^_j.„2(^)- We 
obtain an asymptotic expression of MFPT in terms of the "stochastic po- 
tential" ^{x,V) = — (l/V) In p^y{V). We show in general no continu- 
ous diffusion process can provide asymptotically accurate representations 
for both the MFPT and the PniV) for a DGP. When ni and 112 belong to 
two different basins of attraction, the MFPT yields the T*(y) in terms of 
^{x, V) « (j^oix) + {1/V)(pi{x). For systems with a saddle-node bifurca- 
tion and catastrophe, discontinuous "phase transition" emerges, which can 
be characterized by <I>(x, V) in the limit of y — )■ 00. In terms of time scale 
separation, the relation between deterministic, local nonlinear bifurcations 
and stochastic global phase transition is discussed. The one-dimensional 
theory is a pedagogic first step toward a general theory of DGP. 

1 Introduction 

Nonlinear ordinary differential equations (ODEs) and diffusion processes are two 
important mathematical models, respectively, for dynamics of deterministic and 
stochastic systems. To understand the mathematical properties of these dynamical 
models, it is obligatory to first have a thorough analysis of one-dimensional (1-d) 
systems. In the case of a nonlinear ODE, this is 




(la) 



2 



where x(t) is the state of a system at time t, and in the case of diffusion processes, 

in which u{x,t) is the probability density for a system being in state x at time 
t. A wealth of mathematics has been created by thorough investigations of these 
simple systems. They are now textbook materials with great pedagogic values 
ID 11 El HI. When D{x) = 0, Eq. ^ is reduced to Eq. (dH) via the method 
of characteristics; (flbl) is known as the Liouville equation of (fTal) in phase space. 
The solution to Eq. (flbl) with vanishing D{x) can be considered as a viscosity 
solution to the first-order, hyperbolic partial differential equation. 

In recent years, in connection to mesoscopic size, cellular biochemical dy- 
namics, a new type of mathematical models has emerged: the multi-dimensional 
birth-death process. An A^-dimensional birth-death process is a Markov jump 
process with discrete state i G and continuous time t Q. When applied to 
nonlinear biochemical reaction systems |l6l, its time-dependent probability mass 
distribution, p^t) satisfies the Chemical Master Equation (CME), first studied by 
M. Delbriick, while its stochastic trajectories can be sampled according to the 
Gillespie algorithm [|7l IHI HOl HB . 

The new theory for the Markov dynamics of population systems deserves more 
attentions from applied mathematicians lfT2ll . In addition to its own importance in 
applications, it also provides a unique opportunity for studying the relationship 
between dynamics at mesoscopic and macroscopic levels, which in the past has 
been studied mainly in terms of diffusion processes with Brownian noise. It is a 
widely hold belief that birth-death processes can be approximated by diffusions. 
This turns out not to be the case for nonlinear systems with multiple attractors, as 
we shall show. 

With this backdrop in mind, it is again obligatory to first carry out an com- 
prehensive analysis for a 1-d CME system. Doering et al. have conducted an 
extensive investigation for the asymptotic expressions of the mean first passage 
time (MFPT) |fT3l fT4ll . The aim of the present work is not on this per se, but 
to illustrate the overall mathematical structure of stochastic nonlinear population 
dynamics in terms of the 1-d system. 

The CME for a 1-d birth-death process takes the form 

^Pnit) = Un-lPn-lif) " {Wn + Un)Pn{t) + Wn+lPn+l{t) , ijl > 0) (Ic) 
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in which state-dependent birth and death rates UniV) and WniV) are in general 
functions of n as well as a crucial parameter V, the spatial size or any other ex- 
tensive quantity of the reaction system. For chemical systems consisting of only 
first-order, linear reactions, both n„ and Wn are independent of V ifTSl [T6l [TtII . 
Linear systems have found wide applications in modeling stochastic dynamics of 
single biological macromolecules such as in single-molecule enzymology and 
molecular- motor chemomechanics IfTSl [T9l |20l |2T]| . 

The dependence on V gives rise to a very special feature of the theory of 
Delbriick-Gillespie processes (DGP) and its corresponding CME: One can study 
the important relation between a stochastic dynamical model with a small V and 
a nonlinear deterministic dynamical system with infinitely large V [l22ll23l . T.G. 
Kurtz's theorem precisely establishes such a convergence from the stochastic tra- 
jectories of a DGP to the solution of a nonlinear ODE like Eq. (fTal) . In the 1-d 
DGP, each of (fTal) . (flbl) and (fTcl) has a role. There is also a substantial difference 
between the stochastic system (flbl) in which the stochasticity D{x) and determin- 
istic b{x) are not related per se, while the stochasticity is intrinsic in the dynamics 
of (fTcl) . Therefore models based on diffusion processes are often phenomenologi- 
cal, while the discrete model provides a more faithful representation of a system's 
emerging dynamics based on individual's stochastic behavior. 

Motivated by biochemical applications, recent studies on 1-d DGPs have fo- 
cused on (i) nonlinear bistability, stochastic bistability and multi-stability (i.e., 
multi-modal distribution) lf24H25H26H27l , (ii) non-equilibrium thermodynamics 
and phase transition 112511281 f29ll . (in) large V asymptotics in terms of large de- 
viation theory [l29ll27l . (iv) van't Hoff-Arrhenius analysis motivated by classical 
thermodynamics |l30ll. and (v) dynamics on the circle §^ and oscillations in terms 
of a rotational random walk OTl 13211271 . 

A series of mathematical issues arise in the investigations. The present paper 
initiates a systematic treatment of some of them. One might be surprised by that 
there are still significant unresolved mathematical questions for a one-dimensional 
birth-death process. We simply point out that for large V, the problem under in- 
vestigation is intimately related to the Eq. (ITbl) with a singularly perturbed co- 
efficient D{x) oc (l/V). This is still an active area of research on its own ll33ll . 
In addition, even though straightforward, many explicit formulae in connection 
to the one-dimensional Eq. (ITcl) . also known as hopping models in statistical 
physics, had not been obtained until a need arose from applications. A case in 
point was the 1983 paper of B. Derrida [l34ll . See also 11311141 for recent work on 
the asymptotic analysis of the MFPT problem. Finally, the newly introduced van't 
Hoff-Arrhenius analysis [30] and the analysis for limit cycles [l27l both require 
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Figure 1: Logical schematics showing, for a 1-d DGP, the mathematical relations 
between infinite-time stationary distribution p^*, MFPT for finite-time dynamics, 
and their V oo asymptotics. MFPT for 1-d DGP can be exactly expressed in 
terms of p^* as given in the lower-left box (Eq. [211); Vn ^Iso has an asymptotic form 
shown in the upper-right box. For 1-d continuous diffusion, its stationary density 
function is related to MFPT as shown by the two boxes on the right (Eq. (6]). The 
two MFPTs are "analogous" if we identify Wm with D[^z) and replace summations 
with integrals. A remaining question: What is the asymptotic expression for the 
MFPT in terms of the asymptotic stochasstic potential ^{x). 



consistent asymptotic expansions for large V beyond the usual leading order. 

One of the questions we study in the present work can be succinctly described 
in terms of the diagram in Fig. [T] It is well established that in the limit of large 
V , the stationary solution to the 1-d (flcl) has a WKB (Wentzel-Kramers-Brillouin) 
type asymptotic expansion Pn(V^) ~ exp ( — V <^q{x) — 0i(a;)) where x = n/V 
Il35l [361 |4l. In chemical terms, n is the copy number of a chemical species and 
X is its concentration. Furthermore, it is straightforward to compute the MFPTs 
for both discrete birth-death processes and continuous diffusion [l37ll38l . But it is 
unclear, as indicated by the question mark in Fig. [H whether and how the MFPT 
of the birth-death processes in the limit of large V is related to the "stochastic 
potential function" 0o(x) + {1/V)(f)i{x) obtained from the WKB expansion ll35l 
l36l . This is answered in Eqs. (l25l) and (l34l) . 

For the stationary solution of Kolmogorov forward equation (flcl) . we now have 
a good understanding: For nonlinear dynamical systems with two attractors, there 
is an exponentially large time, e"^ (a > 0) that separates the intra-basin dynamics 
in terms of Gaussian processes IITTll from the inter-basin dynamics of a Markov 
jump process between two discrete states. The "boundary layer" in the singularly 
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perturbed problem is precisely where different attractors join (21 [331. — e~°^ is in 
fact the second largest eigenvalue of the linear system (fTcl) . with zero being the 
largest one. When V" = oo, there is a breakdown of ergodicity |[39ll24ll25l . 

Recognizing this exponentially large time is the key to resolve the so-called 
Keizer's paradox Il24l l40l l4T]| which illustrates the two completely different pic- 
tures for the "steady states" of a deterministic system and its CME counterpart. It 
is also the key to understand the difficulty of approximating a CME like (fTcl ) with 
a diffusion equation like (flbl) for systems with multistability ||39l [25l [TTI . See 
more discussions below. 

The MFPT is the solution to the time-independent backward equation with 
an inhomogeneous term —1 ||37l l38l . An ambiguity arises in the asymptotics of 
MPFT as a WKB solution to the backward equation lfT4ll . This is reminiscent of 
the WKB approach to the stationary forward equation in terms of the nonlinear 
Hamilton-Jacobi equation ll42ll . One of the results in the present work, however, 
is the asymptotic MFPT in relation to the asymptotic stationary solution to the 
corresponding forward equation. 

2 Background on Diffusion Processes 

Because of the intimate relationship between Eqs. (fTcl) and (fTbl) . we shall give a 
brief summary of the relevant results for one-dimensional continuous diffusion in 
Sec. 12.11 Even though it contains no new mathematical result, the presentation is 



novel. Then in Sec. 12.21 we discuss Keizer's paradox from a novel perspective by 
considering a second-order correction to the Kramers-Moyal expansion [|37ll38l . 

2.1 Continuous diffusion 

For a continuous diffusion process with e-small diffusion coefficient: 



plays a central role in its dynamics. In terms of the '^{x), one has the stationary 
distribution 




(2) 



the stochastic potential function 




(3) 



r{x) = Ae 



ivl>(x) 



(4) 
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where A is a normalization factor. is also a Lyapunov function of the ordi- 
nary differential equation dynamics ^ = b{x) since 

Furthermore, the mean first passage time arriving at X2 starting at xi with a re- 
flecting boundary at xq {xq < Xi < X2) is [|37l[38]| 

T.,^.,= Te^^^)/^^ fe-^^yy^dy. (6) 

On the other hand, solving the stationary flux J passing through X2 with Dirichlet 
boundary value f{x2) = 0, leaving the boundary value at xi unspecified, but 
enforcing a normalization condition JJ"^ f^^{x)dx = 1, 



J~' = r e-^^y^^dy r e*(^)/^^. (7) 

Jy D[Z) 

Note that Eqs. ^ and ([7]) are exactly the same if xq = xi in ([61). To understand 
the origin of this intriguing observation, consider the following Gedankenexperi- 
ment: Let a diffusing particle start at xi = xq, which is also a reflecting boundary. 
The particle can only move rightward, and as soon as it hits X2(> Xi), one im- 
mediately takes it back to Xi. Repeating this procedure forms a renewal process. 
Then the mean renewal time is T^^^xj Eq. Now imagine that one connects 
X2 with xi to form a circle, and installs a one-way permeable membrane at the X2- 
Xi junction: a particle that hits from the X2 side goes through the membrane and 
starts at xi instantaneously; but a particle that hits from the xi side is reflected. 
The stationary distribution for the diffusion particle then satisfies f^^{x2) = 0, 
Ixi f^^{'^)dx = 1, and a constant flux J(xi) = J{x2) is the J in Eq. ([7]). 
According to the elementary renewal theorem ^l, Tx^^x2 = J~^- 
Another problem which is widely employed in studies of molecular motor uses 
periodic boundary conditions at xi and X2- Since there is no one-way permeable 
membrane, the boundary condition is f^^{xi) = f^'^{x2) 7^ 0. The cycle flux (i.e., 
mean velocity for a single motor) then is 

J cycle 

The renewal process is then replaced by a semi-Markov process which can go 
both clockwise and counter-clockwise on a circle ||43l . Birth-death processes on 
a circle will be the subject of a forthcoming paper. 
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2.2 Higher-order Kramers-Moyal expansion and Keizer's para- 
dox 



Keizer's paradox was originally introduced to understand a discrepancy between 
the infinitely long time behavior of a CME and its deterministic counterpart in 
terms of nonlinear ordinary differential equations (ODE) Il24l |4T1| . The resolu- 
tion is in the vast separation of time scales: the infinitely long time in the ODE 
is still a very short time in the stochastic dynamics of the CME which involves 
"uphill climbing" and "barrier crossing". The same result also explains the dis- 
crepancy between the stationary distribution of a CME and the stationary distribu- 
tion of the corresponding diffusion approximation via a Fokker-Planck equation 
[|39l |25l . This is now a well-understood subject, intimately related to the finite- 
time condition required in Kurtz's convergence theorem Il22l : The convergence 
when — 7- oo is not uniform with respect to t. 

We now offer a different, more explicit approach to illustrate the diffusion 
approximation problem. The diffusion approximation of a CME in terms of a 
Fokker-Planck equation is actually a truncated Kramers-Moyal expansion up to 
V"^ ||37l[38l . Naturally one can investigate the consequence of keeping the V~'^ 
term in the expansion: 



d 

dx 



d'^u 



du 



e a{x)—-— + eD{x)- b{x)u 



dx 



dx 



0, 



(8) 



where e = 1/V. Applying no-flux boundary condition at x = oo, we have 



c/2 



U 



du 



e^a{x)— + eD{x)- b{x)u = 



(9) 



dx"^ ' ' dx 

We apply the WKB method p4ll2| by assuming the solution to Eq. Q of the form 

1 



u{x) = exp 



X + 



[x] + echoix) + ■ 



(10) 



We then substitute the u{x) into Eq. ^ and collect terms with the leading order 
e° to yield 

a{x) (</.'o(x))' - /^(x)0'o(x) - h{x) = 0. (11) 

We note that if a{x) = 0, then 0o(x) = — j^{h{z) / D{z))dz, as given in Eq. Q. 
In this case, a root of = 0, x*,has </'o(x*) = Oand0Q(x*) = —h'{x*)/D{x*). 
Hence, a stable fixed point of the ODE dx/dt = h{x) corresponds to a local 
minimum of 0o(a^) and a peak in the distribution e~'^°('^)/^ 
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When a{x) ^ 0, Eq. (fTTl) still indicates that at x*, the root of h{x) = 0, 
0q(x*) = 0. Furthermore, differentiating Eq. (fTTl) with respect to x once, we 
obtain 0o(^*) — —b'{x*)/D(x). Therefore, the local behavior of (po{x) near a 
fixed point x* is independent of the higher order terms! However, if a{x) ^ 0, 
then the global behavior of the solution to Eq. (fTTl) will have a non-negligible 
difference from — {b{z) / D{z))dz. This difference contibutes to the difference 

Eq. (fTTl) is in fact a 3rd-order truncated version of the exact equation for 0o(x) 
given by G. Hu[|36l: 

/io(x) (e'^o(^) - l) + Ao(x) (e-<^o(-) - l) = 0, (12) 

with b{x) = ^o{x) — Ao(x), D{x) = (/io(a;) + Ao(x))/2, and a{x) = (Ao(x) — 
fio{x))/Q. The exact, non-trivial, solution to Eq. (fT2l) is dgix) = ln(Ao(x)//io(x)), 
which is given in Eq. (fT6k) below. 



3 One-dimensional Birth-Death Processes: Station- 
ary Distribution and Mean First Passage Time 

We shall now be interested in the Kolmogorov forward equation (fTcl) for the 1- 
dimensional DGP. To be consistent with the macroscopic Law of Mass Action, 
we shall further assume both birth and death rates have asymptotic expansions in 
the limit ofV,n^ oo, n/V — )■ z: 

V-\MV) = /ioW + ^ + ^ + 0(y-^), (13) 
V~'wMV) = Ao(.) + ^ + ^ + +0 (V-^) . (14) 

3.1 Stationary distribution and local behavior near a fixed point 

In terms of the birth and death rates Un{V) and Wn{V), the stationary distribution 
to Eq. (He)) is 
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With large V, applying the lemma in Sec. 17.11 the asymptotic expansion for p^^iV) 
is 

Inp^V(^) - \nf'\x) = V I ln(^44M^ (16) 



Ao(^) 

\i{z) ln(/io(a:)Ao(x)) ^ ^ 



Xo{z) ij,q{z) 

in which we have neglected an x-independent term lnpQ*(V^). 

We shall point out that if one applies a dijfusion approximation to the master 
equation (fTcl) . we obtain the diffusion equation (flbl) with 

D{x) = + (i7a) 



h{x) = /Uo(x) — A 



o(x) + ^(^fii-\i- ^(A'o + fi'o)^ . (17b) 



Then the stationary distribution to the approximated diffusion process is 

Infix) = 2v r ; 7 ; dz m 

Jo mz) + Xo{z) 
^ r Ao(4/ii + X'o - 3fi',) - /xo(4Ai + 3X', + 3fi',) ^ ^ 

Jo (/io + Ao)^ 

which is different from Eq. (fT6l ). even in the leading order. 

However, it is easy to verify that the leading-order terms in the indefinite inte- 
grals in (fT6l) and (fTSi) . as functions of x, have matched locations for their extrema 
as well as their curvatures at each extrema. This is because both 

dlnf'\x) _ y^^ f^ojx) dlnfjx) _ r^y f^oi^) - M^) ^^^^ 
dx Xo{x) dx /Lto(x) + Ao(a;) 

are zero at the root of b{x) = fio{x) — \o{x). One can further check that both have 
identical slopes at their corresponding zeros. 

Therefore, near a stable fixed point a;* of dx/dt = /io(a;) — Ao(x): po{x*) = 
Ao(x*) and p'{x*) < X'{x*), both approaches yield a same Gaussian process with 
diffusion equation 

aj^ ^ ^<^mA _ I _ t)) . (20) 

where ^ = x — x* is widely called fluctuations in statistical physics. This is 
Onsager-Machlup's Gaussian fluctuation theory in the linear regime Il45ll46ll47l 

m. 
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3.2 Mean first passage time and diffusion approximation 



Corresponding to the Kolmogorov forward equation in Eq. (flcl) . the Kolmogorov 
backward equation for the birth-death process is 



dgn 
dt 



WnQn-l - [Un + Wn)gn + ^nfi^n+l- 



(21) 



Then, Tn (0 < n < 71,2), the mean first passage time (MFPT) arriving at n2, 
starting at n with a reflecting boundary at 0, satisfies the inhomogeneous equation 



WnTn-l - {Un + Wn)T„ + M„T„+i = -1, 

with the boundary conditions 



T_i and T, 



n2 



0. 



(22) 



(23) 



The solution can be found in many places, e.g., Ch. XII in ll37l and Equn. 31 
in [[T3l . The result is most compact when expressed in terms of the stationary 
distribution in Eq. (fT5l) : 



T 



n2 m— 1 
m=n+l 1=0 



pf{V) 



(24) 



In the limit of large V, one has the asymptotic expression for T„^_^,j2 (^^^ Sec. 



In 


'm4 

_Mo{z)_ 






Hz) 



V 



in which x = n/V, X2 = n^jV , and 



dz 



In 



My) 



wfa) 

Ao(2/) 



1 



(25) 



^o(x) + -0i(x) + O(r-2^) 



(26) 



given in Eq. (fT6l) . 

Comparing Eqs. (|25l) and Q, we see that the effective "potential function" 



^^{x,V) = ^{x,V) + - In - " ' 



V ln/io(a;) — In Ao(x) 



(27) 
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and effective diffusion coefficient 



D{x) = 

/io(x) 

We note that near fXo{x) = Ao(x), 



In /io(a;) — In Ao(a;) 



(28) 



D{x) ^ Ao(x) 



(Ao - Ho? 



12X1 



(29) 



Following Eq. ©, Eqs. ^ and ^ imply that 

1 - \o{x)/fio{x) 



b{x) 



In yUo(x) — In Ao(a;) 



(/io(x) - Ao(x)) . 



(30) 



As we shall show, even though the general formula for the MFPT given in Eq. (1251) 
has a complex expression for the effective diffusion coefficient D{x) and effective 
drift b(x), the Kramers-like formular for barrier crossing, Eq. (l34l) which only 
involves D(x) at the peak of the potential function, is simple and recognizable. 

Noting the disagreement between the correct asymptotic Eq. (fT6l) and the 
stationary distribution (fTSl ) from the diffusion approximation with D{x) and b{x) 
given in Eq. (flTI) , Hanggi et al. proposed an alternative diffusion equation with 



jjhgt 



fJ'ojx) - Ao(x) 
In Ho{x) — In Ao(a;) 



b{x) = /io(x) - Ao(x), 



(31) 



as a more appropriate approximation for the 1 -dimensional Eq. (fTcl) [39]. While 
the Hanggi-Grabert-Talkner- Thomas diffusion yields the correct leading order sta- 
tionary distribution (fT6l) . we note that the Z}^'^**(x) is different from the D{x) in 
Eq. In fact. 



D''9^\x) In/io -lnAo(x) 



D{x) 
and near fioi^) 



1 - Ao//io 
■ Ao(a;), 

= Ao 



< 1 when — > 1; > 1 when — < 1, 



(32) 



2An 



12X1 



(33) 



The Hanggi-Grabert-Talkner-Thomas diffusion process with diffusion and drift 
given in Eq. (|3T| ) is the only diffusion process that yields the correct deterministic 
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limit b{x) and asymptotically correct stationary distribution (fT6l) . However, the 
discrepancy between D^^^*{x) and D{x) (Eqs. [29]and|33]) indicates that it will 
not, in general, give the asymptotically correct MFPT. Therefore, diffusion pro- 
cesses have difficulties in approximating both the correct stationary distribution 
and the correct MFPT of a birth-death process in the asymptotic limit of large V. 
This conclusion has been dubbed diffusion's dilemma 10111481 . 



3.3 Kramers' formula and MFPT for barrier crossing 

With a correct potential fuunction (I)q{x), the difference between the discrete DGP 
with asymptotic large V and continuous approximation disappears in the compu- 
tations of MFPT for Kramers problem, i.e., barrier crossing. This is because, as 
we have demonstrated, all different approximations can preserve local curvatures 
at the stable and unstable fixed points. At a fixed point x* of h{x): b{x*) = 0, 
and the D{x*) = \o{x*) = fio{x*) in both equations (|29l ) and (l33l) . And with a 
correct "barrier height" and local curvatures at the extrema, the Kramers formula 
is completely determined. 

In fact, applying Laplace's method and considering an energy barrier between 
xl = n\/V and x*2 = n^/V , located at x^, Eq. (l25l) can be simplied into (see Sec. 
Hai for details) 



rji^ rji 



27r 



Mx^W<Pl{x\)\<P'i{xt)\ 

(34) 

in which $(a;,\/) = 0o(a;) + (l/V')0i(x). Note that Eq. OH) contains a (1/V)0i(x) 
term in exponent. This is a key result of the present paper, a new feature for the 
DGP 

Note also that the ambiguity discovered in lfT4]| is associated with MFPT with 
both starting and end points within a same basin of attraction. It does not appear 
in Kramers' formula for inter-basin transition. 

The T* given in Eq. (|34|) plays an all-important role in the dynamics with 
mutiple attractors. It divides the local, intra-basin dynamics from the stochastic, 
inter-basin jump process. The barrier 

V V) - ^xl,V)) = V{Mx^) - 0oK)} + {01 (x*) - 01 (xD), (35) 

thus the time T*, can be increasing or decreasing with V, depending on 0o(x^) — 
00 (x*) > or < 0. This distinction leads to the concept of nonlinear bistability 
vs. stochastic bistability ||26l[T0l[30l . 
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4 Nonlinear Bifurcation and Stochastic Phase Tran- 
sition: A Potential Function Perspective 

Bifurcation is one of the most important characteristics of nonlinear dynamical 
systems [[H |49l. Therefore, a nonlinear stochastic dynamic theory cannot be 
complete without a discussion of its stochastic counterpart. This is still a de- 
veloping area in random dynamical systems |[50l and stochastic processes flU. In 
fact, the very notion of stochastic bifurcation has at least two definitions: the P- 
(phenomenological) and D- (dynamic) bifurcations IISOll . We shall not discuss 
the fundamental issues here; rather provide some observations based on applied 
mathematical intuition. This discussion is more consistent with the P-bifurcation 
advocated by E.C. Zeeman llSTIl . 

The canonical bifurcations in one-dimensional nonlinear dynamical systems 
are transcritical, saddle-node, and pzYc/i/or^ bifurcations [IJ. The saddle-node bi- 
furcation and its corresponding stochastic model have been extensively studied in 
terms of the Maxwell construction ll29l . The key is to realize the separation of 
time scales, and the different orders in taking limits V oo and t — oo. The 
steady states of deterministic nonlinear ODEs are initial value dependent; but the 
steady state distribution of the stochastic counterpart, in the limit of V — >■ oo, is 
unique and independent of the initial value. An ODE finds a "local minimum" of 
^(x) in the infinite time while its stochastic counterpart finds the "global mini- 
mum" at its infinite time. 

In this section. We shall mainly discuss the transcritical bifurcation which has 
not attracted much attention in the past. We shall show that in certain cases, it is 
in fact intimately related to the extinction phenomenon and Keizer's paradox. 

4.1 Transcritical bifurcation 

The normal form of transcritical bifurcation is |[T1 

X = h{x) ^ — X*) — {x — x*Y near x = x* > (36) 

in which the locale of bifurcation is at x = x*; It occurs when /i = 0. Transcritical 
bifurcation is a local phenomenon and the far right-hand-side of (l36l) is the Taylor 
expansion of h{x) in the neighbourhood ofx = x*. Let us assume that the system's 
lower bound is 0; for an ODE to be meaningful to population dynamics, the 6(0) 
has to be non-negative. This implies h{x) has another, stable fixed point x^, < 
x\ < X* when is sufficiently small. 
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Figure 2: Transcritical bifurcation ocurrs between ^2 and at x = x*, as shown 
by the solid and dashed black lines in (A), and corresponding stochastic potential 
shown in (B). Phase transition(s) can occur between fii and fi2, and between 
/is and /i4, when the global minimum of the potential function switches 
between at x* (blue) and near x* (solid black). The global minimun is repre- 
sented by the red curve. The phase transition, therefore, is not really associated 
with the transcritical bifurcation. A very minor "imperfection" can lead the two 
black lines, solid and dashed, to become the two green lines. While the trans- 
critical bifurcation disappeared, the phase transitions are still present. The latter 
phenomenon is structurally stable while the former is not lISTl . The bifurcation is 
local while the phase transitions are global. 

At the critical value of /i = 0, the steady state at x = x* has the form of 
X = —{x — x*y. Hence the corresponding potential function, near x = x*, will 
be 

^.i. - ^ ^0((.-.r). (37) 

Since D{x) > 0, the \l/(x) is neither a minimum nor a maximum at a; = x*. 
There is a minimum of ^(x) at xl. Now for sufficiently small fi ^ 0, a pair of 
minimum and maximum develop in the neighbourhood of x*, approximately at 
X* and X* + /i. Then by continuity, the newly developed minimum of \E'(a;,/i) 
must not he lower than "^{xl, fi). In other words, the stationary distribution of the 
stochastic dynamics, in the limit of V" —t- 00, will not be in the neighbourhood of 
the location of a transcritical bifurcation. 
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However, as in the case of saddle-node bifurcation, phase transition might oc- 
curs for larger Note that the local minimum associated with the transcritical 
bifurcation could become the global minimum of ^{x, fi). Then that occurs, there 
is a phase transition. This is illustrated in Fig. [2l Note however, the phase tran- 
sitions are not associated with the transcritical bifurcation per se: It is really the 
competition between two stable fixed points, the solid green line and the blue line, 

If the X = X* happens at the boundary of the domain of x, we have a more 
interesting scenario. Consider the birth-death system with UniV) = kin and 
WniV) = k_in{n — l)/V + k2n. Then the corresponding iio{x) = kix, Xo{x) = 
A;„ix^ + k2X, and 

b{x) = fio{x) — Ao(x) = {ki — k2)x — x > 0. (38) 

The ODE x = h{x) has a transcritical bifurcation at x = when ki — k2 = 0. 
When ki < /c2, the system has only a stable fixed point at x = 0; when ki > k2, it 
has a stable fixed point at x = \r'^^ , and x = is a unstable fixed point. 

However, the stochastic stationary distribution has a probability 1 at n = 0, 
i.e., extinction, for any value of ki — k2- This is Keizer's paradox ll24l . 

Therefore, when a transcritical bifurcation ocurrs at the boundary of a domain, 
the stochastic steady state exhibit no discontinuous "phase transition". Rather, the 
boundary is an absorbing state. On the other hand, when a transcritical bifurcation 
ocurrs in the interior of the domain, there might not be phase transition associated 
with it. Transcritical bifurcation and phase transition are two different phenomena. 

4.2 Saddle-node bifurcation 

The normal form of saddle-node bifurcation is [[U 

X = h{x) ^ fi — x"^ near x = 0, (39) 

in which the locale of bifurcation is again at x = and it occurs when /i = 
0. Again, it is clear that the stochastic phase transition is not associated with a 
single saddle-node bifurcation event per se. However, it is necessitated by the two 
saddle-node bifurcation events in a catastrophe phenomenon, which has a deeper 
topological root. This is illustrated in Fig. [3l 
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Figure 3: (A) Two saddle-node bifurcations together give the catastrophe phe- 
nomenon. At the four different parameter fi values, the stochastic potential has a 
single minimum (/ii), then passing through saddle-node bifurcation to have two 
minima, with the lower one being the global minimum (^2). Then at n^, the global 
minimum is the upper one. And finally at fi^, saddle-node bifrucation leads again 
to a single steady state. The stochastic phase transition occurs between ^2 and 
Us, denoted by the red dashed line known as the Maxwell construction ll29l . The 
stochastic potential ^(x) corresponding to the four //'s are illustrated below the 
bifurcation diagram. (B) Two saddle-node bifurcations occur at fii and H2. In this 
case, however, there is no phase transition, as illustrated by the given \E'(x) below 
the bifurcation diagram. 

5 van't Hoff-Arrhenius Analysis 

There is a deep connection between the potential function V) and the theory 
of thermodynamics. In this section we provide a brief discussion of the subject, 
which is yet to be fully developed. 

In classical thermodynamics, there is a decomposition of "free energy into 
enthalpy and entropy". The canonical definitions are 




(40) 



(41) 
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It is easy to verfiy that 0o + = "^(^7 Noting an analogue between 

temperature T and 1/^, Eq. POl ) is known as the van't Hoff equation for en- 
thalpy in classic thermodynamics, and Eq. (|4T1) is the definition of entropy: 
S = — {dG/dT)pj^. Then (po and — 0i are the "enthalpic" and "entropic" com- 
ponents of "free energy" $. Knowing one of $(x), (j)o{x), and 0i(x), one can 
determine the other two functions within an additive constant: In classical ther- 
modynamics, they are all called "thermodynamic potentials", and they are equiv- 
alent. 

Since in the limit of V oo, V) = (f)o{x) + + o {V~'^), we 

have _ _ 

lim (I)q{x,V) = (f)o{x), and Um (j)i{x,V) = 0i(x). (42) 

Furthermore, for V^) being a continuous function of V, if 0o(3;, V) = (poi^) 
is independent of V, then l^) = is independent of V. This is because 



(43) 



In this case, is a linear function of 1/K with slope (j)i{x) and intersect 

00 (x). Systems with V independent (f)o(x,V) are "thermodynamically" simple 
systems. 

5.1 Decomposing the potential ^{x,V) into (po{x,V) 

+i^/v)Mx,v) 

The stochastic potential for a one-dimensional birth-death processes has the ex- 
plicit expression: 



xV-l 



Hx,V) = -hnfMV) = ^^\n 

£=0 



Ui{V) - 



(44) 



in which we have neglected, on the right-hand-side, an added term —(1 /V) lnpQ*(V^) 
which is a function of only V , not x. Then following Eqs. (|40|) and (|4TI) . we have 



01 (a;, 



V 



dV 
dV 



+ X 



Vx 



dXnp'^iy) 



dn 



(45) 

n=xV 
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in which we shall substitute the derivative with respect to n by a difference: 



dn 

Thus, the decomposition: 

"n-l 



00 (a;, V) 



d 



^ dV \ 



£=0 
n-l 



Wr. 



(47) 



(48) 

(49) 

n=a;T/ 



5.2 Chemical reaction systems with a linear 1/V dependent 
potential 

What type of CMEs will have a simple, V independent (t)o{x^ V) and V^)? 
Let us consider a simple example: the Poisson distribution with uiiV) = aV, 
Wi{V) = I3£, and 



Pn{V) 



-x*V 



a 



We then have 



xm X + X , 

X* 



(In n\ — n In n + n)^^^y ~ In V2ttx + const. 



(50) 



(51) 

(52) 



Therefore, the Poisson distribution has a stochastic potential with linear l/V de- 
pendence. 

On the other hand, the binomial distribution with u^CV) = k^{etV — i), 
we{V) = k-£, and 



(53) 
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Then, 

0o(a;,1^) 



n-l 



et In 



etV- 

et-x 
et 



X In 



e{etV -n + 1) 



n 



n=xV 



1 Ojet-x) , 
X m h 



n-l ^ , 



X '2V{et — X 

etV -n + 1 



(54) 



n 



+ Inn! + ln(etV^ — n)! 



^=0 " J n=a;y 

fa In a/ (cj — a;)x + const. (55) 

We see that in Eq. (|54|) . there is a (1/V^) dependence. Hence even the simple 
binomial distribution has a ^-dependent (po{x,V). Note that because the po- 
tential function V) defined in Eq. (|44)) is not a continuous function of V, 
V-independent (pi{x, V) no longer dictates V-independent (po{x, V). 



6 Discussion 



Multi-dimensional birth-death processes have become in recent years a funda- 
mental theory for stochastic, nonlinear dynamical systems with individual, "agent- 
based" stochastic nonlinear behavior, and emergent long-time discontinuous stochas- 
tic evolution IfTOl ITTl . The dynamics of such a process is intimately related to 
both nonlinear ordinary differential equations and multi-dimensional diffusion 
processes. In the present paper, we developed a systematic study for the sim- 
plest, one-dimensional system. Each of the equations (fTal) . (flbl) and (fTcl) plays a 
role in the theory of a birth-death process. 

One of the main recent applications of this type of dynamic models is in cellu- 
lar biochemistry in a mesoscopic volume, the Delbriick-Gillespie process (DGP) 
whose Kolmogorov forward equation is widely known as the chemical master 
equation (CME), and whose stochastic trajectories can be sampled following the 
Gillespie algorithm. One special feature of a DGP is a parameter V, the system 
size. When V ^ oc, the trajectory of a DGP becomes the solution to an ODE. 

The application of the CME and the birth-death processes to unimolecular, 
linear reaction systems can be found in lfT5l[T6l[T7l . Steijaert et al. [[52l have also 
presented a coherent summary for the CME with a single variable. In particular, 
they studied the metastability in the thermodynamic limit (V,n ^ oo; n/V = x) 
associated with the Maxwell construction and stochastic phase transition Il28ll29ll . 
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In the present work, we first studied tlie asymptotic stationary solution of the 
1-d stochastic dynamics for large system size V, e~^'^''(a;)-0i(a;) ^^len ex- 
pressed the asymptotic mean first passage time (MFPT) for a 1-d DGP in terms 
of the stochastic potential V) ~ (po{x) + With these two re- 

sults, we obtain an effective diffusion coefficient, Eq. (1281) and a potential. A 
diffusion process defined by these two is significantly different from the diffusion 
approximation proposed by Hanggi et al. [[39ll (comparing Eqs. [28] with [33]) . The 
latter does yield the correct stationary distribution. Together, our analysis shows 
that no single diffusion process with continuous stochastic paths can be globally 
asymptotically accurate for the birth-death processes in general. A situation we 
have called "diffusion's dilemma" |[Tn[48l . 

In the limit of large system size V , a DGP process with multistability (or 
multimodality) exhibits very different dynamics on different time scale: For t -C 
a critical T*, the dynamics is a continuous Gaussian process whose mean value 
follows a deterministic linear dynamic. We call this intra-basin dynamics. For 
t > T*, however, an inter-basin dynamics constitutes Markov jump process that 
moves from an attractor to another attractor. The MFPT gives an estimation for 
the T*, which is in fact the reciprocal of the absolute value of the second largest 
eigenvalue of the system. The largest eigenvalue is always zero for a Markov 
process |[39l[24l[25]l. 

There is a deep relation between the nonlinear bifurcation phenomenon, which 
is when t — )■ oo in a deterministic ODE but still t <^ T*, and stochastic phase 
transition, which is related to t — )■ oo and t ^ T*. In the theory of DGP, phase 
transition can be studied in terms of the stochastic potential V") in the limit 
ofV — 7- oo. The limit of — )■ oo followed by t — oo is widely called quasi- 
steady state which is initial value dependent; while the limit t — )• oo followed 
by A^ — 7- oo is unique, except at the critical point of a phase transition: Phase 
transition occurs in the stationary distribution of an infinitely large system. In the 
present paper for 1-d systems, we have analyzed both transcritical bifurcation and 
saddle-node bifurcation. They are local behaviors while a phase transition is a 
global phenomenon. However, the catastrophe phenomenon necessitates a phase 
transition in terms of the Maxwell construction [[281 (see Fig. [3^). 

Finally, we also suggested an interesting connection between the mathemati- 
cal theory of DGP and the classical thermodynamics in terms of the van't Hoff- 
Arrhenius decomposition of free energy into enthalpy and entropy. 
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7 Mathematical Details 



7.1 A lemma 

Repeatedly using the definition of Riemann integration 

xv-i ^ / f \ 
Jl-Ev^U^)-/ ^(-M-. (56) 

1=0 



we have for a smooth function F{x)\ 

xV-l 











where the coefficients — |, ^, 0, ■ ■ ■ , are the coefficients of the Taylor expansion 



of 

X ^ XX 



2 



i = -± + -^ + ... = yu,x\ (58) 

-1 2 12 ^ ^ ' 

1=1 

Now consider F{x) that is not only a function of x, but also a fuction of V . 
For example ^(^-l)(£-2)/V3 = x{x-l/V){x-2/V) = F{x,V). In the limit 
of y — )■ oo, one has the asymptotic expansion 

F{x, V) = Fo(x) + ^F,{x) + ^F^ix) + ■■■ , (59) 

then 



i=0 



(60) 



7.2 Detailed derivation for Eq. (|25| 

By Taylor expansion and Riemann summation, we shall first show Eq. (l65l) below. 
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We starts with 

m2 m 

F{V) = V' ^ 5^/(r 

m=mi n=0 

s(m+l) /•s{n+l) 

X 



^V[g(m,n)+gx {m,n)Ax+gy{m,n)Ay\ 



dydx, 



(61) 



' s(in) J s(n) 

in which (/a; (x, and Qy (x, denote the partial derivatives {dg / dx) and {dg / dy) 



Furthermore, Ax™ = x — s(m), Ay„ = y — s{n), and s(m) = a + {h — a) 



with (6 — a) = y ^(m2 — mi). Thus 



s{m+l) 



(m) 



(Ax„)^dx 



(Ax, 



+ 1 



s(m+l) 



s(m) 



-£-1 



+ 1 



(62) 



It then can be verified that 

F{V) = ^ ^/(m,n)e^3('"'") 

m=mi 71=0 



E 



/-sCm+l) ;-s(n+l) 



e=0 " Js{m) Js{n) 
m2 m 

m=mi n=0 

°° yl+2 



[gx{m, n)Axm + gy{m, n)Aym) 



dydx 



E 

i=0 



V IV 




>/0 



m2 m 



m=mi n=0 



[gx{m,n)x + gy{m,n)yy dydx 

til 



1+2 




JO 



k=0 



x^g^{m,n)x''gy ^{m,n)y^ '^'^dydx 



m,2 m 



OO £ h / \ f—kt \ 

E E EE^Vw^iTT)! 

m=mi n=0 £=0 k=Q ^ ^ ^ ^ 



m=mi n=0 



g^{m,n)gy{m,n) j^f^^ k\{i-k + 2)\ 
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m2 m 



\n ( 



m=mi n- 



f^gx{m,n)gy{m,n) {i + 2)\ 



X I {gx{m, n) + gy{m, n)Y^'^ - n) - g^^'^{m, n) 



m2 m 



EE 



f{m,n)e^3im,n) 



m=mi 



^ g^{m,n)gy{m,n) I 



o9y{m,n) _|_ ^ 



r?i2 m. 



m=mi n=0 



g,j:{m,n)gy{m,n) 



(63) 



Combining Eqs. (|6TI) and (l63l) we obtain 

Eq. ([63D = = (^^ ^ /(x,?/)e^^(^'^)rf?/dx^ (^V^ ^ 0(V) 



(64) 



Therefore, letting 

/(m,n) = /(m,n) 

we have 

m2 m 



g^{m,n)gy{m,n) 



b px 




a JO 



f{x,y) g, gy e^^("-^ ) 
(gfe - l)(e9!' - 1) 



m=mi n=0 

To show Eq. (|24l) to Eq. (|25l) . we use a = x, 6 = X2, 

■Ao(^) 



dydx I (\/'+0(y)j. 

(65) 



fix, y)e^3ix,y) = g^p 



In 



(iz, and 



Vmx,V)-<l>{y,V)) 
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7.3 Laplace's method for integrals beyond leading order 

One can find the materials in this section from classic texts such as Il44l [2]|. They 
are included here for the convenience of the reader. We first have two important 
formulae: 



2 '•^ 



— / e"* dt = crf{x) 
Jo 



(66) 



and for large x 



Xx-K 



1 ■ 3 ■ 5 ... (2n - 1) 



n=l 



(2x 



2^n 



(67) 



Using Eqs. (1661) and (I67|) . we can evaluate the following integral for z > 0: 



G{x) 



dy 



-(x-z)^/e 



1 - 



7re- 2 



2 , 

"2 22 
(1-2^^) + 



(i-rf.) 



(1-21^) 



2(2;-2)2 



(68) 



X < z 

X = z 

X > z 
(69) 



Near x = z, there is a boundary layer. Let x ^ a = {x — z) / y/e and z ^ rj 
[y — z)/ y/e, then the integral in (|68l) has an inner expression 



G(a) 



_(!t£): 

e 



= -^/e / e ^ dr] 



2 

evr 



erf I — 1 + erf (a) 



;i + erf(a)) 



ee 



2^ 



(2^2 



(70) 
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We have a matched asymptotics: G{0) = G{z), 

Urn Gia) = f (l - = lim G{x) 



2 (T 



2a^ 



lim G{a) 



'ee 



air — 



2a \' 2^2 
Therefore, for a with a unique minimum aty = > 0: 



lim G{x) 



(71) 
(72) 



-In / e 



-V<P{y) 



dy 



-<j){x) + ^ln/;e-^K(")(^-")+^^"(^)(^-^)1rfy X < 
+ iln/;e-^[^<^"("*)(^-"*)'](i?/ X > 

-</>(x)-^ln(\/|0'(x)|)-^^ 



x^l + ^ In — 



Vip" {xi)(^x-xi^ /2 



X < X^ 

X > x"^ 



And within the boundary layer X G (x-'- — \fe,x^ + ^JT) where e = ^2/ (F0"(x^)), 
we have an inner expansion: 



V 



In 



(73) 



With increasing x, the stochastic potential in Eq. (1731) switches from an "enthalpic 
dominant" 

-</'(x)-:^ln(r|0'(x)) + o(V^^), x<x^ 



to a constant independent of x 



— X 



' -' m , . .. , ^, , X > X'''. 



2V V0"(x*) 
The enthalpic term is in fact the limit 



lim - In /" e-^'^^^^rfy = - inf 



(74) 
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Now for the cases in which (j){y) has a single maximum at y = > 0, we 



have: 



^ Jo 



y 



-0(0)-^ In (mO))- 



1 no) 



V i^'(o)f 



< X < 



-0(0) + ^ In 



V(t>'{0) ^ V\<t>'{x)\^ 



-v{cf>{x)-4>{0)) 



X > X 



t 



(75) 



In this case, there is a switching from "entropic dominance" constant when x < 
x^ < X* where (f){x*) = 0(0) to "enthalpic dominance" — 0(a;) when x > x*. 



7.4 Detailed derivation for Eq. (|34| 

For the sake of convenience, we rewrite Eq. (l25l) as follows 



V [ [ f{x,y)e^^^'''yUydx 

J a Jo 



(76) 



where, the same as in Sec. (17.21) . g{x,y) = 0o(a;) — 4>oiy) = ly ^''^j^^^^^- If 
there exists a point {x*,y*) in {{x,y)\xl < x < xl,0 < y < x} that satisfies 
g{x,y) < g{x*,y*),i.e. 



9x{x\y*) = gy{x*,y*) = 0, gxxix*,y*) < 0, gyy{x*,y*) < 0, 
9ly{x*,y*) - gxx{x*,y*)gyy{x*,y*) < 0, 



(77) 



and 



/V,l/*)^0, (78) 
then for large V, denoting Ax = {x — x*) and Ay = {y — y*), 

pb px 

y fix,y)e^3^''^yUydx 

J a Jo 

Ja Jo 

~ yj(a;*^^*)g^^9(a;*,r) r r ^\v[^g,4Axf + ^gyy{Ayf+g^y{Ax)(Ay)]^y^^ 

Jx*—e Jy* — e 
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J — OO J — OO 

/oo /"OO 
-oo J —oo 



-oo J — oo 

27r/(x*,r)e^^("*'^*) ^^^^ 



^Jgxx{x*, y*)gyy{x*, y*) - gly{x*, y*) 
So, from Eqs. (|25l) and (1791) . one can find 



In In ^^4^ 27rp^[<^o(^*'^)-'^o(?'*'^)l+<^i(^*'^)-<^i(2'*'^) 



In In 2Trp^^'f'oi''''y)-My*y)]+M^*y)-My'y) 

Unix*) An U/* I Z,/IC 



27rgy[</.o(x*,y)-<)!.o(y*,v)]+<^i(x*,v^)-<^i(y*,v^) 
27rgF[<i>o(x*,y)-$ofe*,v)] 



) Mo(a;*) y lAto(i/*) Ao(y*) 



Ao(a;*)v/-0o(^*)0o(r)' 



(80) 



where the third equation is because that, at the maximum point {x*, y*), Ao(x*) = 
fioix*), Xo{y*) = /io(y*) and Um^.;.! ^ = 1. The last equation is because Eq. 
implies 0o(x*)0'o'(l/*) < 0. 



References 

[1] L. Perko, Dijferential Equations and Dynamical Systems, 3rd ed., Texts in 
Appl. Math. vol. 7, Springer, New York (2001). 

[2] R. E. O'Malley, Singular Perturbation Methods for Ordinary Dijferential 
Equations, Appl. Math. Sci. vol. 89, Springer- Verlag, New York (1991). 

[3] R. M. Mazo, Brownian Motion: Fluctuations, Dynamics, and Applications, 
Intern. Ser. Monogr. Phys. vol. 112, Oxford Univ. Press, U.K. (2002). 



28 



[4] Z. Schuss, Theory and Applications of Stochastic Processes: An Analytical 
Approach, Appl. Math. Sci. vol. 170, Springer, New York (2010). 

[5] H. Taylor and S. Karlin, An Introduction to Stochastic Modeling, 3rd ed.. 
Academic Press, New York (1998). 

[6] I. R. Epstein and J. A. Pojman, An Introduction to Nonlinear Chemical Dy- 
namics: Oscillations, Waves, Patterns, and Chaos, Oxford Univ. Press, U.K. 
(1998). 

[7] D. A. McQuarrie, Stochastic approach to chemical kinetics, J. Appl. Prob. 
4:413-478(1967). 

[8] D. T. Gillespie, Stochastic simulation of chemical kinetics, Ann. Rev. Phys. 
Chem. 58:35-55 (2007). 

[9] H. Qian and L.M. Bishop, The chemical master equation approach to 
nonequilibrium steady-state of open biochemical systems. Int. J. Mol. Sci. 
11:3472-3500 (2010). 

[10] H. Qian, Cellular biology in terms of stochastic nonlinear biochemical dy- 
namics, J. Stat. Phys. 141:990-1013 (2010). 

[11] H. Qian, Nonlinear stochastic dynamics of mesoscopic homogeneous bio- 
chemical reactions systems - An analytical theory, Nonlinearity 24:R19-R49 
(2011). 

[12] T. G. Kurtz, Approximation of Population Processes, SIAM Pub., Philadel- 
phia, PA (1987). 

[13] C. R. Doering, K. V. Sargsyan and L. M. Sander, Extinction times for birth- 
death processes: exact results, continuum asymptotics, and the failure of 
the Fokker- Planck approximation, Multiscale Modeling Simul. 3:282-299 
(2005). 

[14] C. R. Doering, K. V. Sargsyan, L. M. Sander and E. Vanden-Eijnden, 
Asymptotics of rare events in birth-death processes bypassing the exact so- 
lutions, J. Phys. Cond. Matt. 19:065145 (2007). 

[15] C. Gadgil, C. H. Lee and H. G. Othmer, A stochastic analysis of first-order 
reaction networks. Bull. Math. Biol. 67:901-946 (2005). 



29 



[16] W. J. Heuett and H. Qian, Grand canonical Markov model: A stochas- 
tic theory for open nonequilibrium biochemical networks, J. Chem. Phys., 
124:044110(2006). 

[17] T. Jahnke and W. Huisinga, Solving the chemical master equation for 
monomolecular reaction systems analytically, J. Math. Biol. 54:1-26 (2007). 

[18] H. Qian and E. L. Elson, Single-molecule enzymology: Stochastic 
Michaelis-Menten kinetics. 5?'op/zy5. Chem. 101:565-576(2002). 

[19] H. Qian, Cooperativity and specificity in enzyme kinetics: A single- 
molecule time-based perspective. Biophys. J. 95:10-17 (2008). 

[20] Y. Zhang and M.E. Fisher, Dynamics of the tug-of-war model for cellular 
transport, Phys. Rev. E. 82:011923 (2010). 

[21] Y. Zhang and M. E. Fisher, Measuring the limping of processive motor pro- 
teins, 7. Stat. Phys. 142:1218-1251 (2011). 

[22] T. G. Kurtz, Limit theorems for sequences of jump Markov processes 
approximating ordinary differential processes, J. Appl. Prob. 8:344-356 
(1971). 

[23] J. Keizer, The Mckean model, Kac's factorization theorem, and a simple 
proof of Kurtz's limit theorem. In Probability, Statistical Mechanics, and 
Number Theory: A Volume Dedicated to Mark Kac, G.-C. Rota ed., pp. 1- 
23, Academic Press, New York (1986). 

[24] M. Vellela and H. Qian, A quasistationary analysis of a stochastic chemical 
reaction: Keizer's paradox. Bull. Math. Biol. 69:1121-11 (2007). 

[25] M. Vellela and H. Qian, Stochastic dynamics and nonequilibrium thermo- 
dynamics of a bistable chemical system: The Schlogl model revisited, J. R. 
Soc. Interf. 6:925-940 (2009). 

[26] L. M. Bishop and H. Qian, Stochastic bistability and bifurcation in a 
mesoscopic signaling system with autocatalytic kinase, Biophys. J. 98:1-11 
(2010). 

[27] H. Ge, and H. Qian, Landscapes of non-gradient dynamics without de- 
tailed balance: Stable limit cycles and multiple attractors. Chaos, 22:023140 
(2012). 



30 



[28] H. Ge and H. Qian, Thermodynamic limit of a nonequilibrium steady-state: 
Maxwell-type construction for a bistable biochemical system, Phys. Rev. 
Lett. 103:148103 (2009). 

[29] H. Ge and H. Qian, Nonequilibrium phase transition in mesoscoipic bio- 
chemical systems: From stochastic to nonlinear dynamics and beyond, J. R. 
Soc. Inter/. 8:107-116(2011). 

[30] Y. Zhang, H. Ge and H. Qian, van't Hoff-Arrhenius analysis of transition 
rate dependence upon system's size: Stochastic vs. nonlinear bistabilities in 
population dynamics, h ttp://arXiv.org/abs/ 1 1 1 .2554, (20 1 0) . 

[31] H. Qian, S. Saffarian and E. L. Elson, Concentration fluctuations in a meso- 
scopic oscillating chemical reaction system, Proc. Natl. Acad. Sci. USA 
99:10376-10381 (2002). 

[32] M. Vellela and H. Qian, On Poincar-Hill cycle map of rotational random 
walk: Locating stochastic limit cycle in reversible Schnakenberg model, 
Proc. R. Soc. A. 466:771-788 (2010). 

[33] R. E. O'Malley, Singularly perturbed linear two-point boundary value prob- 
lems, SIAMRev. 50:459-482 (2008). 

[34] B. Derrida, Velocity and diffusion constant of a periodic one-dimensional 
hopping model, J. Stat. Phys. 31:433-450 (1983). 

[35] G. Nicolis and J. W. Turner, Stochastic analysis of a nonequilibrium phase 
transition: Some exact results, PhysicaA 89:326-338 (1977). 

[36] G. Hu, Lyapounov function and stationary probability distributions. Zeit. 
P/jy*. 5 65:103-106(1986). 

[37] N.G. van Kampen, Stochastic Processes in Physics and Chemistry, Rev. and 
enl. ed., Elsevier Science, North-Holland (1992). 

[38] C. W. Gardiner, Handbook of Stochastic Methods for Physics, Chemistry, 
and the Natural Sciences, 2nd ed., p. 263, Springer, New York (1985). 

[39] R Hanggi, H. Grabert, R Talkner and H. Thomas, Bistable systems: Master 
equation versus Fokker-Planck modeling, Phys. Rev. A. 29:371-378 (1984). 



31 



[40] J. Newby and J. P. Keener, An asymptotic analysis of the spatially- 
inhomogeneous velocity-jump process, Multiscale Modeling Simul. 9:735- 
765 (2011). 

[41] P. Childs and J. P. Keener, Slow manifold reduction of a stochastic chemical 
reaction: Exploring Keizer's paradox. Disc. Continu. Dyn. Sys. B 42:1775- 
1794 (2012). 

[42] H. Ge and H. Qian, Analytical mechanics in stochastic dynamics: Most 
probable path, large-deviation rate function and Hamilton-Jacobi equation. 
Int. J. Mod. Phys. B (2012) 

[43] H. Qian and X. S. Xie, Generalized Haldane equation and fluctuation the- 
orem in the steady state cycle kinetics of single enzymes, Phys. Rev. E. 
74:010902(R) (2006). 

[44] J. D. Murray, Asymptotic Analysis, Appl. Math. Sci. vol. 48, Springer- Verlag, 
New York (1984). 

[45] L. Onsager and S. Machlup, Fluctuations and irreversible processes, Phys. 
Rev. 91:1505-1512(1953) 

[46] J. Keizer, On the macroscopic equivalence of descriptions of fluctuations for 
chemical reactions, J. Math. Phys. 18:1316-1321 (1977). 

[47] H. Qian, Mathematical formalism for isothermal linear irreversibility. Proc. 
R. Soc. A Math. Phys. Engr. Sci. 457:1645-1655 (2001). 

[48] D. Zhou and H. Qian, Fixation, transient landscape and diffusion's dilemma 
in stochastic evolutionary game dynamics, Phys. Rev. E 84:031907 (201 1). 

[49] Y. Kuznetsov, Elements of Applied Bifurcation Theory, 3rd ed., Appl. Math. 
Sci. vol. 112, Springer, New York (2004). 

[50] L. Arnold, Random Dynamical Systems, Springer, New York (1998). 

[51] E. C. Zeeman, Stability of dynamical systems, Nonlinearity 1:115-155 
(1988). 

[52] M. N. Steijaert, A. M. L. Liekens, D. Bosnacki, P.A.J. Hilbers and H.M.M. 
ten Eikelder, Single-variable reaction systems: Deterministic and stochastic 
models, Ma?/z. Biosci. 227:105-116 (2010). 



32 



